function [ flowX flowY ] = HornSchunk(imageSeq, params)
    params.preSmoothingWindow = 2 * ceil(norminv(0.9 , 0, params.preSmoothingSigma));
    params.numScales = 2 + floor(log(min(size(imageSeq{1}, 1), size(imageSeq{1}, 2)) / 16));

    for iIm = 1 : numel(imageSeq) - 1
        im1 = im2double(MyRGB2Gray(imageSeq{iIm}));
        im2 = im2double(MyRGB2Gray(imageSeq{iIm + 1}));         
                 
         pyramid1 = ImagePyramid(im1, params);
         pyramid2 = ImagePyramid(im2, params);
         
         u = zeros(size(pyramid1{params.numScales}));
         v = u;
         for iScale = params.numScales : -1 : 1             
             [ u v ] = ProcessScale(pyramid1{iScale}, pyramid2{iScale}, u, v, params);
             
             if (iScale > 1)
                 u = imresize(u, size(pyramid1{iScale - 1}));
                 v = imresize(v, size(pyramid1{iScale - 1}));           
             end
             
             if params.isVerbose
                fprintf(1, 'Horn Schunk: Image %d/%d | Scale %d/%d Done.\n', iIm, ...
                    numel(imageSeq) - 1, params.numScales - iScale + 1, params.numScales);
             end
         end         
        
        if (numel(imageSeq) > 2)        
            flowX{iIm} = u;
            flowY{iIm} = v;
        else
            flowX = u;
            flowY = v;
        end
        
        if params.isVerbose
            fprintf(1, 'Horn Schunk: Image %d/%d Done.\n', iIm, numel(imageSeq) - 1);
        end
    end
end

function pyramid = ImagePyramid(image, params)
    pyramid = { image };
    h = fspecial('gaussian', params.preSmoothingWindow, params.preSmoothingSigma);    
    for i = 2 : params.numScales
        curImage = imfilter(pyramid{i - 1}, h, 'symmetric');
        curImage = imresize(curImage, 0.5);
        pyramid{i} = curImage;
    end
end

function [ u v ] = ProcessScale(im1, im2, u, v, params)    
    for iWarp = 1 : params.numWarpIterations         
        [ x y ] = meshgrid(1 : size(im1, 2), 1 : size(im2, 1));            
        x2 = x + u;
        y2 = y + v;
        im2Warped = interp2(x, y, im2, x2, y2);        
        
        flow = [ u(:) ; v(:) ];

        [ lhs rhs ] = HornSchunkEquations(im1, im2Warped, flow, params);
        dFlow = lhs \ rhs;
        flow = flow + dFlow;            

        u = reshape(flow(1 : numel(flow) / 2), size(im2));
        v = reshape(flow(1 + numel(flow) / 2 : end), size(im2));        

        if (params.medianFilterSize > 0)
            u = medfilt2(u, params.medianFilterSize);
            v = medfilt2(v, params.medianFilterSize);                       
        end
    end
end